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This  paper  deals  with  the  development  of  guidance,  navigation  and  control  algorithms  for  a 
prototype  of  a  miniature  aerial  delivery  system  capable  of  high-precision  maneuvering  and  high 
touchdown  accuracy.  High  accuracy  enables  use  in  precision  troop  resupply,  sensor  placement,  urban 
warfare  reconnaissance,  and  other  similar  operations.  Specifically,  this  paper  addresses  the  terminal 
phase,  where  uncertainties  in  winds  cause  most  of  the  problems.  The  paper  develops  a  six  degree-of- 
freedom  model  to  adequately  address  dynamics  and  kinematics  of  the  prototype  delivery  system  and 
then  reduces  it  to  a  two  degrees-of-freedom  model  to  develop  a  model  predictive  control  algorithm  for 
reference  trajectory  tracking  during  all  stages.  Reference  trajectories  are  developed  in  the  inertial 
coordinate  frame  associated  with  the  target.  The  reference  trajectory  during  terminal  guidance,  just 
prior  to  impact,  is  especially  important  to  the  final  accuracy  of  the  system.  This  paper  explores  an 
approach  for  generating  reference  trajectories  based  on  the  inverse  dynamics  in  the  virtual  domain. 
The  method  results  in  efficient  solution  of  a  two-point  boundary-value  problem  onboard  the  aerial 
delivery  system  allowing  the  trajectory  to  be  generated  at  a  high  rate,  mitigating  effects  of  the 
unknown  winds.  This  paper  provides  derivation  of  the  guidance  and  control  algorithms  and  present 
analysis  through  simulation. 


I.  Introduction 

MANEUVERABLE  ram-air  parafoils  are  widely  used  today.  The  list  of  users  includes  skydivers,  smoke 
jumpers,  and  Special  Forces  to  name  a  few.  Furthermore,  their  extended  range  compared  to  round  parachutes 
makes  them  very  practical  for  payload  delivery.  For  round  canopy  parachutes  the  current  requirement  for  a  drop 
zone  (DZ)  size  for  the  High  Velocity  Container  Delivery  System  from  2,500m  altitude  as  defined  in  the  Joint  U.S. 
Army  /  U.S  Air  Force  Field  Manual  5-430-00-21  is  as  large  as  780m  by  1,610m  or  a  1,255,800m2  (260  acres) 
footprint.  During  the  latest  humanitarian  aid  aerial  delivery  missions,  such  as  the  delivery  of  food,  supplies  and 
medications  to  war-torn  areas  like  Bosnia,  Sarajevo,  Afghanistan,  and  Iraq  50%  of  payloads  were  never  recovered. 
Such  requirements  for  parafoils  haven’t  been  developed  because  high-glide  ratio  delivery  systems  cannot  glide  to  a 
specific  location  by  themselves,  but  rather  require  a  guidance,  navigation  and  control  (GNC)  unit  to  produce  and 
track  the  corresponding  steering  commands.2  3 

Recent  introduction  of  the  Global  Positioning  System  (GPS)  made  the  development  of  fully  autonomous  ram-air 
parafoils  possible.  Moreover,  ram-air  parafoils  are  considered  a  vital  element  of  the  unmanned  air  vehicles 
capability,4"6  can  be  used  in  conjunction  with  a  space  station  crew  rescue,7"10  and  in  some  other  applications." 
Autonomous  parafoil  capability  implies  delivering  the  system  to  a  desired  landing  point  from  an  arbitrary  (limited 
by  winds,  deployment  condition  and  system  performance)  release  point  using  an  onboard  computer,  sensors  and 
actuators.  The  navigation  subsystem  manages  data  acquisition,  processes  sensor  data,  and  provides  guidance  and 
control  subsystems  with  information  about  parafoil  states.  Using  this  information  along  with  local  wind  profiles,  the 
guidance  subsystem  plans  the  mission  and  generates  a  feasible  trajectory  to  the  desired  landing  point.  Finally,  it  is 
the  responsibility  of  the  control  system  to  track  this  trajectory  using  the  information  provided  by  the  navigation 
subsystem  and  onboard  actuators. 

The  need  for  development  of  the  fully  autonomous  different-weight  aerodynamic  decelerator  systems  emerged 
just  in  the  past  decade.  Natick  Soldier  Center,  Natick,  MA  pursues  and  pushes  the  development  of  the  new 
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generation  joint  precision  aerodynamic  delivery  systems  (JPADS)  capable  of  achieving  a  breakthrough  accuracy  of 
less  the  100m  circular  error  probable  (CEP).12  The  underlying  concept  to  be  employed  on  every  system  is  to  have  a 
field-laptop-based  JPADS  Mission  Planner  (developed  by  Planning  Systems,  NOA/FSL  and  DL)  that  can  predict  the 
winds  at  different  altitudes  with  unprecedented  accuracy.13"15  Wind  prediction  enables  computation  of  the  required 
release  point  from  the  airdrop  cargo  airplane.  Improved  Container  Delivery  System  developed  for  26’  Ring  Slot  or 
G-12  unguided  ballistic  parachutes  along  with  JPADS-Mission  Planner  allows  reduction  of  the  DZ  to  109m  by 
219m  or  23,937m2  (5  acres),  i.e.  52  times  smaller  then  for  the  conventional  systems  that  do  not  have  the  JPADS 
Mission  Planner.  The  authors  were  lucky  to  participate  in  the  program  of  the  development  of  the  Affordable  Guided 
Aerodelivery  Systems16  that  not  only  employed  JPADS-Mission  Planner,  but  also  by  relatively  inexpensive  means 
allowed  conversion  of  regular  ballistic  (dumb)  circular  parachutes  into  a  guided  system. 17  The  accuracy  achieved 
with  this  system  is  even  smaller  than  5-acre-footprint  and  remains  on  the  order  of  potentially  more  accurate  parafoil- 
based  systems.18 

The  ultimate  goal  for  the  perspective  delivery  system  is  to  have  payloads  delivered  from  the  large  standoff 
distances  with  <100m  CEP  accuracy.  So  during  the  last  decade,  several  GNC  concepts  for  gliding  parachute 
applications  have  been  developed  and  published.19'31  Specifically,  the  best  systems  were  demonstrated  during 
Precision  Airdrop  Technology  Conference  and  Demonstration  PATCAD-2001,  PATCAD-2003,  PATCAD-2005 
and  PATCAD-2007  at  U.S.  Army  Yuma  Proving  Ground,  Yuma,  AZ,  and  Precision  Airdrop  Demonstration 
Capability  PCAD-2006  and  PCAD-2008  near  Toulouse,  France.32"34 

Analyzing  the  results  of  dozens  of  airdrops  it  can  be  stated  that  the  most  accurate  systems  can  only  achieve 
100m-150m  accuracy  (excluding  SPADES  by  Dutch  Space25  that  assures  about  40m-60m  accuracy).  More  detailed 
analysis  reveals  that  GNC  systems  developed  by  Charles  Stark  Draper  Laboratory,”  "  Georgia  Tech”  "  and  others 
implement  simple  proportional-integral-derivative  (PID)  controllers  to  track  the  reference  heading,  which  is 
assigned  based  on  a  combination  of  some  logics  and  heuristic  rules.  For  the  final  stage  of  descent  to  the  DZ  some 
controllers  use  a  moth  mode,  others  attempted  to  utilize  a  precomputed  landing  trajectories  data  base.24  Though 
some  of  the  systems  utilize  JPADS  Mission  Planner  they  only  use  it  to  compute  a  release  point,  but  they  do  not  use 
its  winds  prediction  during  the  following  control. 

This  paper  builds  on  the  generic  GNC  algorithm  developed  earlier27'28  and  specifically  concentrates  on  the 
terminal  phase  of  guided  descent,  where  the  majority  of  errors  contribute  to  the  final  touchdown  accuracy.  The  paper 
is  organized  as  follows.  Section  II  overviews  the  overall  control  strategy  consisting  of  six  phases  that  follow  each 
other  upon  loosing  altitude  during  a  descent.  Section  III  features  the  Model  Predictive  Control  (MPC)  used  to  track 
the  desired  trajectory  at  each  phase.27  28  Section  IV  introduces  a  six  degree-of-freedom  model  of  a  generic  parafoil 
delivery  system  and  characterizes  the  linearized  model,  developed  for  use  with  MPC.  Section  V  provides  details  of 
the  simple  straightforward  control  algorithm  during  the  final  phase  (turning  base),  while  Section  VI  deals  with  the 
more  sophisticated  adaptive  optimal  control  for  this  final  phase.  The  paper  ends  with  several  computer  simulations 
and  conclusions. 


II.  Guidance  Strategy 

In  the  current  design  five  parameters  need  to  be  defined  by  the  user  before  deployment.  The  five  parameters  are: 
target  location,  away  distance,  cycle  distance,  turn  diameter,  and  wind  heading  angle.  These  five  parameters  define 
five  fixed  tracking  points,  which  are  the  Target,  A,  B,  C,  and  D  as  shown  in  Fig.l.  Using  these  five  points,  the  ADS 
precision  placement  objectives  are  defined.  As  shown  in  Fig.l,  trajectory  planning  for  precision  placement  is 
divided  into  six  phases.  They  are  as  follows: 

Phase  0:  After  canopy  opening,  guidance  is  delayed  to  allow  oscillations  to  cease; 

Phase  1:  The  initial  phase  is  used  to  reach  the  loitering  area  (bounded  by  points  A,  B,  C,  and  D),  where  it 
is  assumed  that  the  system  is  released  upwind  of  the  loiter  area; 

Phase  2:  Loitering  is  commenced.  The  system  estimates,  winds,  descent  rate,  and  the  altitude  at  which  to 
turn  towards  the  target.  The  loiter  area  is  defined  by  four  points  A,  B,  C,  and  D.  As  stated  above 
these  four  points  are  determined  by  the  target  location,  wind  direction,  away  distance,  cycle 
distance,  and  turn  diameter.  During  loiter,  winds  are  estimated  when  traveling  from  A  to  B  and 
from  C  to  D  (long  parallel  legs).  Ideally,  at  least  one  complete  loop  around  A,  B,  C,  and  D  is 
completed.  Winds,  distance  to  the  target,  and  vertical  velocity  are  used  to  determine  a  switching 
altitude  zstart  to  start  the  downwind  leg  towards  the  target  (as  explained  in  Section  V).  When  the 
measured  altitude  reaches  zstan  Phase  2  is  terminated  and  the  system  turns  toward  the  target; 

Phase  3:  Wind  estimation  is  halted  as  the  loiter  area  is  exited; 
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Phase  4:  Wind  estimate  is  reinitiated  as  the  system  travels  toward  the  target.  All  estimated  data  is  used  to 
determine  a  distance  Dswitch  past  the  target  to  make  a  final  turn  and  approach.  When  the  distance 
Dswitch  is  reached  Phase  5  is  initiated; 

Phase  5:  A  180-degree  turn  to  the  final  approach  glide  slope; 

Phase  6:  Final  approach  and  landing  into  the  wind. 


Fig.l  Guidance  strategy. 

The  precision  placement  algorithm  tracks  a  desired  heading  using  a  model  predictive  control  (MPC)  described  in 
Sections  III  and  IV.  In  each  phase  the  current  position  and  desired  tracking  points  are  used  to  determine  a  desired 
heading  of  the  system.  Prior  to  terminal  guidance  (phases  0  through  4)  the  desired  heading  is  found  simply  by 
heading  directly  to  the  next  tracking  point.  However,  in  the  critical  final  phases  the  desired  heading  is  found  by 
rapidly  calculating  an  optimal  trajectory,  detailed  in  Sections  V  and  VI. 

III.  Model  Predictive  Control 

Consider  a  simple  Single-Input  Single-Output  (SISO)  discrete  system  described  in  state  space  form  as: 

x*+i  =Axa-+b«a-  (|) 

yk  =cxk 

(Hereinafter  bold  font  for  lower-case  symbols  denotes  vectors  and  bold  font  for  upper-case  symbols  denotes 
matrices).  Assume  that  the  system  matrices  A,  B,  and  C  are  known  and  that  xk  is  the  state  vector,  uk  is  the  control 

input,  and  yk  is  the  output  at  time  k .  The  model  described  above  can  be  used  to  estimate  the  future  state  of  the 
system.  Assuming  a  desired  trajectory  wk  is  known  an  estimated  error  signal  ek  =  wk  —yk  is  computed  over  a  finite 
set  of  future  time  instants  called  the  prediction  horizon,  H  (hereinafter  the  symbol  “A”  is  used  to  represent  an 

estimated  quantity).  In  model  predictive  control,  the  control  computation  problem  is  cast  as  a  finite  time  discrete 
optimal  control  problem.  To  compute  the  control  input  at  a  given  time  instant,  a  quadratic  cost  function  is  minimized 
through  the  selection  of  the  control  history  over  the  control  horizon.  The  cost  function  can  be  written  as: 

y  =  (W-Y)rQ(W-Y)  +  UrRU,  (2) 

where 

w  =[wk+uwk+2,...,wk+Hp]r ,  (3) 

Y  =  KcA+Kc^U,  (4) 

U  =  [uk,uk+\,...,uk+H'_t  1  ,  (5) 

and  both  R  and  Q  are  symmetric  positive  semi-definite  matrices  of  size  H p  x  II p  .  Equation  1  is  used  to  express  the 
predicted  output  vector  Y  in  terms  of  the  system  matrices 
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CA 

CA2 

K  Cil  =  :  ,  (6) 

_CAHp  _ 

CB  0  0  0  0  " 

CAB  CB  0  0  0 

Kcab  =  CA2B  CAB  CB  0  0.  (7) 

:  :  :  0 
CA^P  B  CA2B  CAB  CB 

The  first  term  in  Eq.(2)  penalizes  tracking  error,  while  the  second  term  penalizes  control  action.  Equations  (2), 
(6),  and  (7)  can  be  combined  resulting  in  the  cost  function 

T  =  (W-KcA-Kc^U)rQ(W-KcA-Kc^U)  +  UrRU  (8) 

that  is  now  expressed  in  terms  of  the  system  state  \k  ,  desired  trajectory  W,  control  vector  U  and  system  matrices  A, 
B,  C,  Q,  and  R. 

The  control  U,  which  minimizes  Eq.(8),  can  be  found  analytically  as 

U  =  K(W-KcA),  (9) 

where 

K  =  {^cabQ^-cab  +  R)  *  KqtsQ  •  (10) 

Equation  10  contains  the  optimal  control  inputs  over  the  entire  control  horizon,  however  at  time  k  only  the  first 
element  uk  is  needed.  The  first  element  uk  can  be  extracted  from  Eq.(9)  by  defining  K,  as  the  first  row  of  K.  The 
optimal  control  over  the  next  time  sample  becomes 

«*=K1(W-KC4it),  (11) 

where  calculation  of  the  first  element  of  the  optimal  control  sequence  requires  the  desired  trajectory  W  over  the 
prediction  horizon  and  the  current  state  \k  . 

IV,  Parafoil  Modeling 

The  full  combined  flexible  system  of  the  parafoil  canopy  and  the  payload  can  be  represented  as  a  9  or  8  degree- 
of-freedom  (DoF)  model  depending  on  the  specific  harness  connecting  these  two  pieces  together.  Alternatively  it 
can  also  be  modeled  as  a  solid  structure  requiring  only  six-DoF,  which  include  three  inertial  position  components  of 
the  system  mass  center  as  well  as  the  three  Euler  orientation  angles  of  the  parafoil-payload  system.  This  section  first 
introduces  a  six-DoF  model  of  a  generic  parafoil  system  that  will  be  used  for  simulation  and  GNC  algorithm 
development,  and  then  addresses  a  simplified  model  used  to  compute  the  MPC  gains  appearing  in  Eq.(l  1). 

A.  Six  DoF  Model 

Figure  2  shows  a  schematic  of  a  parafoil  and  payload  system.  With  the  exception  of  movable  parafoil  brakes,  the 
parafoil  canopy  is  considered  to  be  a  fixed  shape  once  it  has  completely  inflated.  A  body  frame  {,5}  is  fixed  at  the 
system  mass  center  with  the  unit  vectors  iB  and  k/;  orientated  forward  and  down.  Orientation  of  the  parafoil 
canopy  with  respect  to  the  body  frame  is  defined  as  the  incidence  angle  F. 

Orientation  of  the  body  frame  {B}  is  obtained  by  a  sequence  of  three  body-fixed  rotations.  Starting  from  the 
inertial  frame  {/},  the  system  is  successively  rotated  through  Euler  yaw  y/,  pitch  9,  and  roll  <p.  The  resulting 
transformation  from  the  inertial  to  body  frame  is 

cect//  cesy  ~se 

T/B  =  S0S0Ci/s —C0Sy/  S0S0' Sy/  +  C0Cl//  S</>C0  >  (12) 

c</>se( y  +  s0st//  c0s0sp  -  s,/>ci//  c<j>c0  _ 
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where  the  common  shorthand  notation  for  trigonometric  functions  is  employed,  so  that  for  any  angle  a  sin(a)  =  su 
and  cos(a)  =  ca  .  Translation  kinematics  of  the  system  become 


(13) 


where  x,  y,  and  z  are  the  inertial  positions  and  u,  v,  and  w  are  body  frame  velocities.  Defining  the  angular  velocity  of 
the  parafoil  payload  system  as: 

™B/I  =  PiB+Cl}B+rkB  (14) 

results  in  the  following  rotational  kinematics  equations: 


X 

u 

y 

II 

.H 

&  *-5 

V 

z 

w 

0 

1  sAe  cAd 

P 

e 

- 

0 

1 

q 

(15) 

V 

0  s<!>ce  c<t>ce 

r 

( ta  =  tan(a) ). 


Fig.2  A  generic  6-D0F  parafoil-payload  system. 

Forces  and  moments  acting  on  the  parafoil  and  payload  have  contributions  from:  weight,  aerodynamic  loads,  and 
apparent  mass  of  the  canopy.  Weight  contribution  of  the  system  is  expressed  in  the  body  frame. 


F,= 


~se 

s<t>ce 

C^Cn 


mg , 


(16) 


with  m  being  the  system  mass  and  g  the  acceleration  from  gravity.  Aerodynamic  forces  and  moments  have 
contributions  from  both  the  canopy  and  payload.  Both  canopy  and  payload  contributions  are  combined  into  a  single 
aerodynamic  model  using  standard  aerodynamic  derivatives.  The  aerodynamic  force  is  modeled  as: 

2 " 


F  =  T 


pv2asp 


AB  ' 


a 


Q)0  +  C Da  2 

Cypfi 

Q.O  +  C Laa 


(17) 


where;  a  and  /?  are  the  angle  of  attack  and  sideslip  of  the  body  frame,  Va  is  the  true  airspeed,  Sp  is  the  parafoil 
canopy  area,  and  T/)w  is  the  transformation  from  the  aerodynamic  to  body  frame  by  the  angle  of  attack  a.  Similarly, 
the  aerodynamic  moment  is  written  as: 
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M„  = 


pvX 


biCv+^C^  +  ^r  +  d-'q 


ISa 


8 a) 


c  ( CmO  +  Cmaa  +  4“  ^ ma ?) 


*(C, 


nP  +  2Vn 


2  V  mq'i 

b  Cnp P+  2y  Cnrr  +  d'CnSaSa) 


(18) 


where,  b  and  c  are  the  canopy  span  and  chord,  Sa  is  the  asymmetric  brake  deflection,  and  d  is  the  maximum 
brake  deflection.  Parafoil  canopies  with  small  mass  to  volume  ratios  can  experience  large  forces  and  moments  from 
accelerating  fluid  called  “apparent  mass.”  They  appear  as  additional  mass  and  inertia  values  in  the  final  equations  of 
motion.  Parafoil  canopies  with  small  arch-to-span  ratios  and  negligible  camber  can  be  approximated  to  useful 
accuracy  by  an  ellipsoid  having  three  planes  of  symmetry;  however,  the  planes  of  symmetry  in  the  canopy  frame 
may  not  be  aligned  with  the  body  frame,  as  shown  in  Fig.2.  Apparent  mass  forces  and  moments  for  an 
approximately  ellipsoidal  canopy  can  be  conveniently  defined  using  the  translation  and  angular  velocities  expressed 
in  the  canopy  frame  \P\.  Defining  T BP  as  the  single  axis  transformation  from  the  body  to  canopy  reference  frame 

by  the  incidence  angle  r,  rBM  =  [xBM  ,yBM ,zBM^  as  the  vector  from  the  system  mass  center  to  the  apparent  mass 

center,  and  W  as  the  wind  vector,  the  velocity  of  the  canopy  at  the  apparent  mass  center  can  be  expressed  in  the 
canopy  frame  as: 


=  T„ 


/r u 1 


w 

\LVVJ 


+sf 


XBM 

y  bm 

ZBM 


-T«W 


(19) 


where  the  second  term  represents  the  vector  cross  product  using  the  skew-symmetric  matrix  Sf ,  constructed  as 
follows: 


S'4  - 
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0 
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(20) 


(here  4*,  4  an(J  4  are  the  components  of  vector  4  expressed  in  a  coordinate  frame  {A}). 
Similarly,  the  angular  velocity  expressed  in  the  canopy  frame  {P}  becomes 


p 

P 

q 

r 

II 

q 

r 

(21) 


Forces  and  moments  from  apparent  mass  and  inertia  are  then  found  by  relating  the  fluid’s  kinetic  energy  to 
resultant  forces  and  moments.35  Apparent  mass  and  inertia  contributions  expressed  in  the  body  frame  can  be  written 
as: 


Ffl.m.  -  T bp 
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where 


l 

WJ 

~  A 

0 

0" 

~P 

o 

o 

^a.m. 

0 

B 

0 

and  I  a.i.  = 

0 

Q  o 

0 

0 

c 

0 

0  R 

(24) 


(Apparent  mass  and  inertia  values  A,  B,  C,  P,  Q,  and  R  appearing  in  Eq.(24)  can  be  calculated  for  known  simple 
shapes  or  can  be  approximated  as  discussed  in  Ref.35.) 
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Dynamic  equations  are  formed  by  summing  forces  and  moments  about  the  system  mass  center,  both  in  the  body 
reference  frame,  and  equating  to  the  time  derivative  of  linear  and  angular  momentum.  Final  dynamic  equations  of 
motion  are  expressed  compactly  in  matrix  form  below,  where  for  the  quantities  in  Eq.(24)  the  common  convention  is 
used  for  tensors  of  second  rank  such  that  I \  =  TjPI  JTbp  : 


u 

V 

wI3x3  +l'a.m. 
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i  +  r,  -sf  i'  sf 

p 

M 

rBM  a-m ■ 

a-L  rBM  a-m •  rBM  _ 
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where  I3x3  is  the  identity  matrix,  I  = 
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(25) 


(26) 


(27) 


B.  Linearized  Model 


The  guidance  strategy  outlined  in  Section  II  requires  a  simple  SISO  controller  with  yaw  as  an  output  and  the 
brake  deflection^  as  the  control.  A  simple  2-DoF  model  of  the  roll  and  yaw  dynamics  is  used  since  for  parafoils 
pitch  and  speed  are  not  typically  controllable.  The  state  vector  x  for  the  2-DoF  rotational  model  contains  roll,  yaw, 
roll  rate,  and  yaw  rate: 


x  =  [<j>,i//,p,r]T .  (28) 

Equations  (15)  and  (25)  describe  the  nonlinear  rotational  kinematics  and  dynamics.  However,  for  MPC  (1)  a 
linear  discrete  model  is  required.  Assuming  that  the  aerodynamic  velocity  Va  is  constant,  Eqs.(15)  and  (25)  can  be 
numerically  linearized  about  the  steady  state 


"  o" 
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(29) 


With  the  convention  s*  =s-s0  t  the  resulting  linear  discrete  model  for  the  parafoil  with  parameters  listed  in  the 
Appendix  and  sampling  period  of  0.5s  is: 
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Matrices  given  by  Eq.(30),  along  with  the  output  matrix  C  =  [0  1  0  0]  define  the  optimal  brake  deflection 
according  to  the  MPC  algorithm  in  Eq.(l  1). 
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T  - 1 

±turn  *1 


Tapp  ~  h  h 


V.  Terminal  Guidance 

The  last  three  phases  introduced  in  Section  II  (Phases  4-6)  are  the  most  critical  stages  of  parafoil  guidance.  The 
system  must  be  dropped  up  wind  of  the  target  to  ensure  it  can  be  reached,  however,  it  is  desired  to  impact  the  target 
traveling  into  the  wind  to  reduce  ground  speed.  In  addition  it  is  beneficial  to  arrive  near  the  target  with  excess 
altitude  in  order  to  make  final  guidance  decisions.  Finally,  there  is  a  very  strict  time  limitation.  The  ADS  can  be 
slightly  late  or  earlier  departing/arriving  to  all  other  phases,  but  this  last  Phase  6  ends  sharply  at  landing.  All  this 
means  that  special  precautions  have  to  be  made  in  building  a  control  algorithm  for  the  terminal  phase. 

First,  assume  that  there  is  no  cross  wind  component,  i.e.  that  the  ground  winds  uploaded  to  the  system  before 
deployment  have  not  changed.  An  ideal  terminal  guidance  trajectory  can  be  defined  as  outlined  in  Fig. 3,  where  the 
following  notations  are  used  (Fig. 3  represents  the  left  approach  pattern,  but  everything  will  be  the  same  for  the  right 
pattern  as  well.): 

x,y,z  -  a  standard  right-handed  North-East-Down  coordinate  frame  with  the  origin  at  the  target; 

tstart  -  time  corresponding  to  the  beginning  of  the  downwind  leg  (Phase  4); 

t0  -  time  corresponding  to  the  beginning  of  the  final  180°-turn  (Phase  5); 

tj  -  time  corresponding  to  intercepting  the  final  approach  (Phase  6); 

texjt  -  time  when  guidance  switches  from  Phase  5  to  Phase  6; 

t2  ~  time  of  touchdown; 

L  -  distance  away  from  the  target  line; 

^switch  ~  optimal  distance  to  pass  the  target  (initiating  a  final  turn  at  t0  should  achieve  the  desired 

impact  location); 

-t0  -  final  turn  time; 

-  final  approach  time,  determining  t0  and 

T) switch  (a  large  value  allows  correction 
for  terminal  errors  and  reduces  errors 

from  changing  winds,  while  a  small  value  _ 

reduces  errors  during  final  approach);  t 

R  -  final  turn  radius;  /J  , 

W  -  wind  (positive  when  coming  from  the  Xj  i 

direction  and  negative  as  shown  in  Fig. 2);  t 

(//([)  -  final  turn  function  for  parafoil  to  track  £ 

during  the  final  turn;  [ 

D  =  WTturn  -  distance  defined  by  imperfection 
(asymmetry)  of  the  final  turn  because  of 
the  wind. 

The  terminal  guidance  problem  can  be  summarized  as 
follows.  For  a  parafoil  in  the  presence  of  wind  W,  at  altitude  z  , 

and  a  distance  L  from  the  target,  find  the  distance  Dswitch  to  the  final  turn  initiation  point  (TIP),  required  to  travel 
past  the  target  for  an  ideal  impact  at  t2  (of  course,  Dswitch  is  a  function  of  the  desired  .) 

In  general,  the  dynamic  model  of  a  parafoil  is  complex  and  nonlinear,  so  that  the  problem  can  only  be  solved 
numerically.  However,  in  what  follows  we  will  make  some  assumptions,  allowing  an  analytical  solution  to  be  used 
as  a  reference  trajectory  in  the  control  algorithm.  These  assumptions  are:  a)  the  turn  rate  is  slow,  so  that  the  roll  and 
sideslip  angles  can  be  ignored,  and  b)  the  descent  rate  V*  and  air  speed  F;*  are  viewed  as  nearly  constant  (defined 
by  the  canopy,  lines,  etc.  and  treated  as  the  known  values).  In  this  case  it  immediately  follows  that 

h=t start  -“pr-  •  (31) 

The  problem  reduces  to  handling  a  simple  kinematic  model  represented  by  three  components  of  the  ground 
velocity  as: 


Fig.3  Terminal  guidance  maneuver. 
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X 

-W  +  V,*  cos  y/ 

y 

= 

Vl  sin  i// 

z 

V* 

r  V 

(32) 


Now,  the  reference  turn  function  ys(t)  can  be  chosen  as  any  function  that  satisfies  the  boundary  conditions 
y/{to)  =  0  and  y/[tx )  =  -n  (for  the  left  pattern  as  shown  on  Fig. 3).  For  instance,  for  a  linear  turn  (constant  turn  rate) 
the  solution  to  the  required  turn  rate  and  the  final  turn  time  Tturn  are  provided  by: 


V, 


71 R 


¥=-—  and  Tturn=  —  . 


R 


V, 


(33) 


Flence,  the  most  straightforward  algorithm  to  control  the  descending  system  at  the  terminal  phase  is  to  control  its 
turn  rate  for  example  as  follows: 


Vh 

v  =~+T 


(34) 


(the  plus-minus  signs  correspond  to  the  left  and  right  turn,  respectively).  After  the  final  turn  the  system  travels 
directly  to  the  desired  target.  Assuming  that  the  wind  W  is  constant  all  way  down  and  a  constant  turn  rate  (33), 
integration  of  inertial  velocities  along  axes  x  and  y  from  tstarl  to  t2  (Phases  4-6)  yields  two  simple  equalities: 

to 


D. 


switch 


-  j*  xdt  -  Ji 


-  xdt  =  D, 


switch 


-WTturn-(Vh  +  W)T  =  0 , 


(35) 


z  +  K, 


L  +  D 


switch 


Vh  -  W 


+  V  t  +V  t  =0 

1  r  v  ‘'turn  1  r  v  lapp  w  ' 


Resolving  them  with  respect  to  Dswitch  and  a  T  results  in 


^ 'switch  WI t 'turn 


fv;2-w2) 

( 

— -T 

L+wrtum\ 

{  2V,  J 

*  turn 

V  v 

vh-w  ) 

T  = 

app 


Vh -W 
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—  -T 

T/*  turn 
J\Vv  J 


L  +  WThl 

2v: 


(36) 

(37) 

(38) 

and  Tn„„  become.  As  the 


app 


From  Eqs.(37)  and  (38)  it  can  be  seen  that  the  higher  the  altitude  z  ,  the  larger  Dswjtch 
parafoil  loiters  upwind  of  the  target,  zstart  the  altitude  at  which  to  turn  towards  the  target  can  be  found  by  using  a 
desired  final  approach  time  T%Z  .  The  switching  altitude  to  achieve  is  then  given  by  solving  Eq.(38)  for  z 
leading  to  the  following  expression: 


z  =  -V 

*  start  '  v 


T  L  +  WThim 

,um  v;-w 


2Vh 


v„-w 


■pdes 

Lapp 


des 


=  V, 


L  +  Vh  ( Ttum  +  2Tapp 


) 


w-v. 


(39) 


V  h  h  J 

Once  the  system  is  traveling  towards  the  target  the  goal  is  to  bring  the  system  from  its  initial  state  at  the  top  of 
the  downwind  led  to  the  point  defined  by  xT  =  Dswitch  and  yT  =  +2R  (for  the  left  and  right  turn,  respectively).  The 

distance  Dswitch  is  estimated  during  the  downwind  leg  constantly  using  the  analogue  of  Eq.(37): 


D switch  ^VTtum  + 


f  v;2-w2) 

f-z  T  x  +  WTturn] 

z{Vl2  -W2)+  v;vxm (Vl  -W)  +  xv; (Vl  +  W) 

l  2Vh  J 

V*  turn  y*w 

V  v  r  h  ,r  J 

2  KK 

,(40) 


where  Vh  ,  Vv  and  W  are  the  estimates  of  the  corresponding  parameters  at  current  position  (x,z  ).  Note,  Eq.(40) 
produces  the  value  of  Dswitch  in  the  assumption  that  V*h  ,  V*  and  W  remain  constant  from  the  current  altitude  all 
way  down. 

Figure  4  demonstrates  simulation  results  for  the  ideal  case  when  all  parameters  (horizontal  airspeed,  descent  rate, 
wind)  are  assumed  to  be  estimated  with  no  errors,  there  are  no  disturbances,  and  the  commanded  yaw  rate  is  tracked 
precisely.  Specifically,  for  horizontal  airspeed  of  22.4ft/s,  descent  rate  of  lOft/s,  wind  of -lOft/s,  turn  radius  of  125ft, 
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and  the  desired  final  approach  time  of  7.5s,  the  maneuver  starts  at  the  altitude  of  561.36ft  (Z=750ft).  The 

distance  Dswitch  is  estimated  to  be  -82.3ft,  it  takes  17.5s  to  make  a  full  turn,  and  the  system  touches  down  precisely 
on  target  in  45.65s. 


Fig.4  Simulated  guidance  in  ideal  conditions:  3D  trajectory  (a),  horizontal  projection  and  the  time  history  of  the  yaw 

angle  (b). 


In  practice,  sensor  errors,  variable  winds  (magnitude  and  direction),  and  imperfect  control  will  certainly  disturb 
the  ideal  touchdown  depicted  above.  For  example.  Table  1  provides  measurement  errors  typical  for  the  GNC  unit 
used.  Winds  can  change  (and  are  changing)  all  the  time.  Let  us  illustrate  how  inaccuracies  in  measuring  ADS 
position  errors  arriving  at  the  TIP  and  variable  winds  may  affect  the  touchdown  accuracy. 


Table  1.  Errors  in  measuring  the  ADS  states. 


Parameter 

Bias 

Standard  deviation 

Altitude  z 

10  ft 

1  ft 

Descent  rate  z 

0  ft/s 

2  ft/s 

Roll  angle  (j) 

1.5  0 

1.5  0 

Y aw  angle  y/ 

3  ° 

1.5  0 

Roll  rate  p 

1  Ys 

1.5  °/s 

Yaw  rate  r 

1  °/s 

1.5  °/s 

Figure  5  shows  Monte  Carlo  simulations  where  two  cases  were  considered.  First,  it  was  assumed  that  the 
downwind  leg  ends  up  with  some  error  in  the  ADS  horizontal  position  and  heading.  The  standard  deviation  of 
horizontal  position  error  in  each  coordinate  was  assumed  to  be  20ft  and  the  standard  deviation  in  heading  10  °.  As 
seen  from  Fig. 5a  these  errors  result  in  spreading  the  arrival  to  the  final  approach  point.  The  situation  is  even  worse 
when  unaccounted  winds  are  acting  on  the  ADS  during  the  final  turn. 

For  simulations  shown  in  Fig. 5b,  the  standard  deviation  for  wind  disturbances  applied  at  each  integration  step  of 
0.05s  was  only  3ft/s  (with  a  zero  mean  value)  in  both  x  and  y  direction,  but  as  shown  it  led  to  a  drastic  degradation 
of  the  overall  performance.  Obviously,  if  we  add  a  constant  (unmeasured  wind)  that  acts  on  the  system  all  way  down 
in  any  direction  it  will  simply  blow  the  ADS  off  the  desired  trajectory.  This  is  the  worst  case  scenario.  Now  the 
question  is  how  to  mitigate  the  effect  of  changing  winds? 

The  control  algorithm  only  provides  indirect  estimates  of  the  wind  based  on  the  difference  between  the  parafoil’s 
measured  ground  speed  and  its  steady  state  speed  components,  V*  and  V*h  .  Again,  these  estimates  are  only  valid  for 
the  current  altitude  and  provide  no  prediction  of  future  winds  at  lower  altitudes.  Therefore,  rather  than  trying  to 
better  estimate  current  winds,  the  authors  proceeded  with  the  development  of  a  robust  algorithm  accounting  for 
constantly  varying  winds  during  the  final  turn  and  compensate  for  them  by  changing  the  yaw  rate  command.  The 
next  section  addresses  this  algorithm  based  on  the  inverse  dynamics  in  the  virtual  domain  in  more  detail. 
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Fig.5  Final  turn  trajectories  with  not  precise  initial  conditions  (a)  and  unaccounted  winds  (b). 


VI.  Optimal  Terminal  Guidance 

To  overcome  the  difficulties  of  fighting  unaccounted  winds  the  following  two-point  boundary-value  problem 
(TPBVP)  was  formulated.  Staring  at  some  initial  point  at  t  =  0  with  the  state  vector  defined  as  x0  =  [x0, yQ,Yof 

we  need  to  bring  our  ADS  influenced  by  the  last  known  constant  wind  w  =  [lf,0, 0]r  to  another  point. 


X/  = 


\vl:+w) 


t£,o,-x 


(for  consistency  with  the  previous  section  we  will  consider  a  left  pattern  turn)  at 


t  =  tj  .  Hence,  we  need  to  find  the  trajectory  that  satisfies  these  boundary  conditions  along  with  a  constraint  imposed 
on  control  (yaw  rate),  \y/\  <  (/>max  ,  while  finishing  the  maneuver  in  exactly  Tturn  seconds  (if  needed,  the  performance 
index  may  also  include  additional  terms  related  to  the  compactness  of  the  maneuver  and/or  control  expenditure).  The 
found  optimal  control  y/op,(t)  is  then  tracked  by  the  ADS  GNC  unit  using  the  MPC  algorithm  in  Sections  III  and  IV. 

Unaccounted  winds  wrf,s/  =  \wx,  wv,oJ  (the  vertical  component  will  be  accounted  for  indirectly  via  estimating  the 

current  Vv )  will  not  allow  exact  tracking  of  the  calculated  optimal  trajectory.  Computation  of  the  optimal  trajectory 
will  be  updated  during  the  final  turn,  each  time  starting  from  the  current  (off  the  original  trajectory)  initial  conditions, 


requiring  the  ADS  to  be  at  x.-  = 


\vl  +  w) 


t£,o,-x 


in  Tturn  seconds  from  the  beginning  of  the  turn.  The 


remaining  time  until  arrival  at  the  final  approach  is  computed  as 


tf  =  -^r~TK 


final 


(41) 


where  Tfif  is  the  final  calculated  approach  time  from  Eq.(38)  before  entering  the  final  turn. 

Development  of  the  optimal  trajectory  algorithm  begins  by  recalling  from  Eq.(32)  the  horizontal  trajectory 
knimematics: 


X 

_y_ 

(42) 


-W  +  Vh  cos  y/ 

V*h  sin  y/ 

Recognizing  that  if  the  final  turn  horizontal  trajectory  is  given  by  Eq.(42)  the  yaw  angle  along  this  trajectory  is 
related  to  the  change  of  inertial  coordinates  as  follows: 

-i  y 


y/  =  tan 


(43) 


x  +  W 

Differentiating  (43)  provides  the  yaw  rate  control  required  to  follow  the  reference  final  turn  trajectory  in  presence  of 
constant  wind  W: 
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(44) 


¥  = 


y(x  +  W)-xy 


(x  +  Wf+f  ' 

The  inertial  (ground)  speed  along  the  trajectory  will  also  depend  on  the  current  yaw  angle: 

\VG |  =  y/x2  +  y2  =  Vf;2  +W2- 2V*Wcosi/s  .  (45) 

Now,  following  the  general  idea  of  direct  methods  of  calculus  of  variations  (e.g.,  see  Ref.36)  we  will  assume  the 
solution  of  the  TPBVP  to  be  represented  analytically  as  functions  of  some  scaled  abstract  argument 
r  =  r/zy  e  [0;  1] : 

x(f )  =  Pi(t)  =  al0  +  a\r  +a\f 2  +  a\r2  +b\  sin {xf)  +  b\  sm(2xr), 
y(r )  =  P2(t)  =  a0  +a[r  +a2z~  +a3r  +b{  sin(^r)  +b2  sin(2^r). 

The  coefficients  a?  and  ( r]  =  1, 2  )  in  these  formulas  are  defined  by  the  boundary  conditions  up  to  the  second- 
order  derivative  at  r  =  0  and  r  =  zy  (  f  =  1 ).  According  to  the  problem  formulation  and  Eq.(42)  these  boundary 
conditions  are  as  follows: 


(46) 
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_y_ 
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0 

(47) 

(48) 


While  the  final  conditions  (48)  will  be  constant  (the  second  derivatives  are  zeroed  for  a  smooth  arrival),  the 
initial  conditions  will  reflect  the  current  state  of  the  system  at  each  cycle  of  optimization.  Let  us  now  differentiate 
Eqs.(46)  two  times  with  respect  to  r  to  get: 

TfP'{r )  =  a'l  +2alz +3a1r2  +nb1  cos(^r)  +  27rbj  cos(2;zt), 

(49) 

r2P”(j)  =  2al  +6a^f  - tt2^  sin(^zr)-(2^-)2&2  sin(2 jtf). 

Equating  these  derivatives  at  the  terminal  points  to  the  known  boundary  conditions  (47)-(48)  yields  a  system  of 
linear  algebraic  equations  to  solve  for  coefficients  a(  and  (r/  =  1, 2  ).  For  instance,  for  the  x-coordinate  we’ll  get 


(50) 
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Being  resolved  this  system  yields 


=  x0,  a\=-(x0-xf)  + 


(2xq  +x"f)rj 


a1 


a,=- 


(xq  ~x"f)r2f 


U1  _  2(x'o -x'f)Tf+(x o  +x"f)r2f 
b\  ~ 


j  12(x0  —Xf)  +  6(xq  +x'j-)tj-  +(xq  -x"f)r2 


(51) 


b\  = 


4/r  24^- 

The  only  problem  is  that  the  derivatives  in  Eqs.(51)  are  taken  in  the  virtual  domain,  while  the  actual  boundary 
conditions  are  given  in  the  physical  domain.  Mapping  between  the  virtual  domain  [  0;  t  f  \  and  physical  domain 

[0 \tf  \  is  addressed  by  introducing  a  speed  factor  X: 


i  dr 
4-  -  —  > 
dt 


(52) 


Using  this  speed  factor  we  may  now  compute  corresponding  derivatives  in  the  virtual  domain  using  the  obvious 
differentiation  rules  valid  for  any  time -variant  parameter  <f: 
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Inverting  Eqs.(53)  yield: 


(53) 


£  =  *4',  £  =  *(*'?+ A4'). 


?  =  r'£,  %"  =  r2£-A,'rl£ .  (54) 

Note  that  we  only  need  to  use  Eqs.(54)  once  to  transfer  the  boundary  conditions.  Since  the  speed  factor  X  simply 
scales  the  entire  problem  -  the  higher  speed  factor  A  is,  the  larger  zy  it  results  in  (see  explanation  in  Ref.  3  7)  -  we 

may  let 

A0.j  =1  and  =  0  ,  (55) 

which  means  we  can  safely  assume 

?  =  £,?  =  £.  (56) 

Now  let  us  describe  the  numerical  procedure  for  finding  the  optimal  solution  among  all  candidate  trajectories 
described  by  Eqs.(46).  First,  we  guess  the  value  of  the  only  varied  parameter  zy  and  compute  the  coefficients  of  the 

candidate  trajectory  using  Eqs.(51)  with  the  boundary  conditions  (47)  and  (48)  converted  to  the  virtual  domain  via 
Eqs.(54)  (accounting  for  Eqs.(55)).  For  the  initial  value  of  r  r  we  can  take  the  length  of  the  circumference 

connecting  terminal  points: 

Tf  +07- Fo)2  •  (57) 

Having  an  analytical  representation  of  the  candidate  trajectory,  Eqs.(46)  and  (49),  define  the  values  of  x  - , 
yj,x'j,  and  y'j,  j  =  \,...,N  over  a  fixed  set  of ./V points  spaced  evenly  along  the  virtual  arc  [0;zy]  with  the  interval 

AT  =  Tf(N-iy\  (58) 

so  that 

Tj=Tj_l+AT,  j  =  2,...,N  ,  (n  =  0  ).  (59) 

Then,  for  each  node  j  =  2,...,  N  we  compute 


A0- 1=. 


(xf-x,- 1)2 +(y,-y,- \f 


'V*2  +W2  -2Vj*Wcosy/,  l 


(  ¥\  =  ¥0  X  and 

Aj  =  ArA  . 

The  yaw  angle  <//  can  now  be  computed  using  the  virtual  domain  version  of  Eq.(43): 

-1  ^y'j 


y/  i  =  tan 


Ajx'j  +  W 


(60) 

(61) 

(62) 


Finally,  the  yaw  rate  \j/  is  evaluated  using  Eq.(44)  (with  time  derivatives  evaluated  using  Eq.(53))  or  simply  as 

yj=(y/j-y/j_x)Atf_l.  (63) 

When  all  parameters  (states  and  controls)  are  computed  in  each  of  the  N  points,  we  can  compute  the  performance 


index 


where 


J  = 


N- 1 


V2 


I>y4 


+  k^A , 


A  =  max(0;|^.|-^.max)2 


(64) 


(65) 


with  k v  being  a  scaling  (weighting)  coefficient.  Now  the  problem  can  be  solved  say  in  the  Mathworks’ 
development  environment  with  as  simple  built-in  optimization  function  as  fminbnd.  It  should  be  noted  that  the 
numerical  algorithm  based  on  the  direct  method  introduced  in  this  section  with  the  non-gradient  optimization  routine 
fminbnd  based  on  the  straight  forward  golden  section  search  and  parabolic  interpolation  algorithm  allows  computing 
the  optimal  turn  trajectory  very  fast.  To  be  more  specific  a  16bit  80MHz  processor,  allowed  computation  of  the  17.5 
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second  turn  maneuver  with  N  =  20  in  10  iterations,  which  took  only  0.07s  all  together.  With  the  controls  update  rate 
of  0.25s  that  means  that  the  trajectory  can  be  updated  as  often  as  every  control  cycle. 

As  a  demonstration  of  the  proposed  approach,  Fig. 6  shows  an  example  of  the  optimal  final  turn  trajectory,  which 
is  essentially  the  same  as  in  Fig.4b.  Flowever,  as  opposed  to  Fig. 4b  the  yaw  angle  changes  smoothly  at  departure 
from  TIP  and  arrival  to  the  top  of  the  final  approach.  As  a  result  the  turn  trajectory  is  slightly  different  from  that  of 
Fig. 4b,  but  what  is  of  the  most  importance  is  that  ADS  captures  the  final  approach  glide  slope  in  exactly  Tturn 
seconds  as  Fig. 4b  trajectory  does. 

Therefore,  now  we  have  a  tool  allowing  us  to  construct  the  optimal  trajectory  from  any  initial  point  to  the 
predetermined  final  point  to  be  achieved  in  a  certain  time.  To  this  end,  even  if  we  get  to  the  TIP  with  errors  as 
shown  in  Fig. 4a,  we  can  easily  accommodate  these  errors  and  provide  the  optimal  control  so  that  we  will  still  be  at 
the  top  of  the  final  approach  in  the  predetermined  time,  Tturn  (see  Fig. 7). 


Fig.6  Optimal  guidance  in  the  perfect  conditions  as  in  Fig.4b.  Fig.7  Optimal  guidance  with  an  error  at  TIP. 


Similarly,  in  case  of  wind  disturbances  (Fig.5b)  the  ability  to  recompute  the  trajectory  from  the  current  (off  the 
original  trajectory)  conditions  to  the  same  final  conditions  as  often  as  practical  seems  to  alleviate  the  problem.  To 
this  end,  Fig.8  presents  two  examples  of  handling  wind  disturbances  with  the  trajectory  update  rate  set  as  a  function 
of  tracking  error.  Once  the  tracking  error  while  following  the  previously  generated  trajectory  exceeds  a  certain 
amount  (10ft  and  15ft  in  Fig. 8a  and  8b,  respectfully)  a  new  trajectory  is  being  generated.  In  Fig.8  multiple  reference 
trajectories  are  shown  with  the  solid  lines,  while  a  continuous  parafoil  trajectory  is  shown  with  a  dashed  line.  The 
points  where  the  old  reference  trajectory  was  abandoned  are  depicted  with  the  circles,  and  the  corresponding 
locations  of  the  parafoil  blown  away  from  the  reference  trajectory  by  strong  (up  to  15%  of  the  forward  speed) 
unaccounted  winds  is  marked  with  triangles.  In  practice,  the  accuracy  of  a  GPS  sensor  is  on  the  order  of  10ft,  so 
rather  than  making  trajectory  updates  according  to  a  sliding  error  as  shown  in  Fig.8,  it  is  also  possible  to  time 
schedule  updates. 

Although  the  optimal  terminal  guidance  algorithm  works  well,  several  further  adjustments  need  to  be  made  in 
order  to  make  it  more  robust  in  practice.  The  kinematic  models  in  Eq.(42)  does  not  account  for  parafoil  turning 
dynamics  and  assumes  that  the  sideslip  and  roll  angles  are  small.  When  the  turn  rate  is  sufficiently  small  or  the 
radius  R  is  large,  the  model  in  Eq.(42)  provides  sufficient  accuracy.  Flowever,  in  order  to  track  the  desired  trajectory 
for  a  wide  range  of  R,  the  error  from  sideslip  and  turning  dynamics  can  be  compensated  by  adding  an  additional 
commanded  i//  for  the  first  tpre  seconds  after  the  TIP  has  been  reached.  Specifically,  for  the  first  tpre  seconds  of  the 
maneuver  the  command  yaw  rate  y/‘  produced  by  Eq.(62)  is  augmented  as  follows: 


,~c  _  v 

V  V  K-turn  D 

K 

where  Klvrn  is  the  correction  gain  and  Vh  j R  is  the  required  constant  turn  rate  in  (33). 


(66) 


Disturbances  from  wind  and  sensor  errors  during  tracking  of  the  optimal  trajectory  will  result  in  errors  in  final 
approach  arrival.  These  errors  can  be  accommodated  by  re -computing  the  optimal  trajectory  from  the  current  (off 
the  original  trajectory)  conditions  to  the  desired  final  approach  point  during  Phase  5.  Re-computation  of  the  optimal 
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trajectory  is  scheduled  so  that  the  specified  number  of  updates  occurs  at  equal  intervals  over  the  period  from  t0  to 

texit' 


Fig.8  Optimal  guidance  with  wind  disturbance  with  a  10ft  (a)  and  15ft  (b)  cap  on  a  tracking  error. 

The  final  capture  location  Xf(x  at  tf  in  Eq.(48))  assumes  that  during  phase  6  the  parafoil  travels  directly  to  the 
target  by  the  shortest  path.  In  practice  disturbances  result  in  deviation  from  the  shortest  path  and  the  guidance 
algorithm  will  result  in  a  trajectory  with  less  than  perfect  efficiency  at  reaching  the  target.  In  order  to  correct  for 
realistic  errors  during  the  final  approach  the  final  capture  location  xyis  corrected  to  be 

(67) 

where  £j-maI  is  the  final  approach  efficiency. 

Figures  9  and  10  show  examples  of  how  these  additions  to  the  guidance  law  work  for  desired  turn  radii  of  125ft 
and  325ft,  respectively.  In  Figs.  9  and  10  the  full  nonlinear  parafoil  model  is  simulated  with  the  MPC  control 
algorithm  described  previously  used  to  track  the  desired  y/  angle.  In  both  cases  Tapp  is  7.5s,  £rlnai  is  0.95,  tt-texi,  is 

3.0s,  tpre  is  6.0s,  Ktmn  is  1.0,  Vh  is  22.4ft/s,  Xf  is  88.4ft,  and  the  wind  W  is  -lOft/s.  Results  are  shown  from  the  TIP 
(x0,  y0,  z0)  for  both  cases  which  are  (-88.3ft,  250ft,  -31 1ft)  for  Fig.9  and  (-363ft,  650ft,  -658ft)  for  Fig.  10.  As  seen, 
even  with  a  single  trajectory  update  (Fig. 9a)  the  parafoil  lands  near  desired  location  while  additional  updates  allow 
for  more  correction  near  the  end  of  the  trajectory  (Fig. 9b).  In  Fig.  10  the  longer  turn  allows  more  error  to  accumulate 
with  even  one  trajectory  update  resulting  in  a  sufficient  correction. 

VII.  Precision  Placement  Results 

The  complete  precision  placement  algorithm  is  developed  by  combining  the  loitering  pattern  (Phases  1  and  2), 
exit  decision  (Phase  3)  based  on  terminal  guidance  geometry,  estimation  of  TIP  (Phase  4),  and  the  optimal  terminal 
guidance  and  the  final  approach  (Phases  5  and  6).  The  parafoil  loiters  in  Phase  2  continuously  estimating  wind  W  , 
air  speed  Vh  ,  and  descent  rate  V* .  The  altitude  at  which  to  exit  loitering  and  start  towards  the  target  zstart ,  found  by 
using  Eq.(39)  assumes  that  the  parafoil  is  already  traveling  towards  the  target.  While  loitering  in  Phase  2  there  will 
be  a  delay  between  reaching  the  altitude  zstart  and  when  the  parafoil  can  turn  to  face  the  target.  In  order  achieve  the 

desired  final  approach  time  T^pp  ,  the  parafoil  must  exit  loitering  Tdel  seconds  early.  The  altitude  to  exit  loitering 
can  then  be  estimated  by  using  Eq.(39)  using  the  combined  time  Tdel  +  T^pp  rather  than  the  ideal  T^pp  .  Upon 
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exiting  the  loiter  area  Phase  4  is  entered  and  the  parafoil  travels  toward  the  target  updating  estimates  of  Dswitch  .  The 
location  where  the  parafoil  reaches  Dswitch  is  defined  as  the  TIP  and  the  beginning  of  Phase  5.  After  exiting  the  loiter 
pattern  the  approach  time  T  is  determined  by  Eq.(36).  Due  to  disturbances,  measurement  error,  and  tracking  error 

T  may  not  match  T  .  The  optimal  final  trajectory  is  determined  using  the  last  estimates  of  W  ,  Vh  ,  Vv  and 
T  as  outline  in  Section  IV.  While  in  Phase  5  the  optimal  final  turn  trajectory  is  continually  updated  until  time  texit , 
after  which  Phase  6  begins  the  final  approach. 


Fig.9  Optimal  real-dynamics-augmented  guidance  with  one  (a)  and  three  (b)  trajectory  updates  with  a  125ft  turn  radius. 


Fig.10  Optimal  real-dynamics-augmented  guidance  with  one  (a)  and  three  (b)  trajectory  updates  with  a  325ft  turn  radius. 


The  full  algorithm  is  demonstrated  through  numerical  integration  of  the  nonlinear  model  starting  from  (x ,y,z)  of 
(-1200ft,  -200ft,  -2200ft).  All  precision  placement  parameters  used  are  provided  in  Table  2.  Winds  are  1  Oft/s 
downrange  (W  =  -lOft/s)  both  at  altitude  and  ground  level.  During  the  final  turn  in  Phase  5  two  updates  of  the 
optimal  trajectory  are  used.  In  all  phases  of  guidance  the  MPC  algorithm  based  on  the  linear  parafoil  approximation 
is  used  to  track  the  desired  y/  angle. 

Results  of  the  algorithm  are  shown  in  Figs.  11  to  13.  Progression  of  algorithm  phases  is  shown  by  triangles 
representing  the  transition  to  the  next  guidance  phase,  with  the  first  transition  marking  the  beginning  of  Phase  2.  The 
trajectory  is  shown  in  Fig.  1 1  with  the  four  loiter  region  points  marked  by  x’s,  optimal  trajectory  updates  marked  by 
squares,  and  the  final  target  by  a  circle.  Upon  entering  Phase  2  the  parafoil  circles  the  loiter  region  in  a  clockwise 
direction. 
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Table  2.  Parameters  of  the  precision  placement  algorithm. 

Parameter 

Value 

Away  distance 

750  ft 

Cycle  distance 

400  ft 

Final  turn  radius,  R 

125  ft 

rpdes 

1  app 

7.5  s 

T 

delay 

6.5  s 

400 

V  j  max 

20  °/s 

tpre 

6.0  s 

1 1  ~  t exit 

3.0  s 

8 final 

0.95 

K-turn 

1.0 

Figure  12  shows  the  estimation  of  the  exit  altitude  zstart  from  Eq.(39)  based  on  current  wind  estimates,  the  approach 
time  T^s  ,  and  transition  delay  Tdel  .  Transition  from  Phase  2  to  3  occurs  when  the  altitude  reaches  zstart  at  120 
seconds  at  an  altitude  of  732ft.  Figure  13  shows  that  Tapp  is  14.0s  ( T‘‘pp  +Tdelay  )  when  Phase  3  is  entered.  Transition 
to  Phase  4  occurs  at  133s  at  which  time  Tapp  has  decreased  to  6.9  seconds  due  to  the  time  required  to  make  the  turn 
exiting  the  loitering  region.  During  Phase  4  it  can  be  seen  that  the  required  Tapp  continues  to  vary  as  tracking  error 
and  wind  estimates  change.  In  Phase  4  the  transition  to  the  final  turn  occurs  when  Dswitch  ,  the  estimated  distance  to 
pass  the  target,  is  reached.  Figure  12  shows  ~Dswitch  ,  where  the  value  is  positive  when  the  parafoil  must  turn  before 
reaching  the  target  and  negative  if  the  turn  is  to  occur  after  passing  the  target.  Early  in  the  trajectory  -Dswitch  is 
much  less  than  zero  demonstrating  if  a  turn  to  target  was  performed  at  that  time  the  parafoil  would  need  to  pass  the 
target  by  a  large  amount  before  turning  because  of  it  current  altitude.  Transition  from  Phase  4  to  5  occurs  when  the 
distance  to  target  L  reaches  -Dswitch  at  154  seconds  with  the  parafoil  62ft  ahead  of  the  target.  Two  updates  of  the 
optimal  final  turn  trajectory  occur  before  transition  to  Phase  6  at  169  seconds.  Final  impact  occurs  at  2ft  cross  range 
and  8ft  down  range  at  180.0  seconds. 
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Monte  Carlo  simulations  of  50  drops  were  completed  using  the  complete  precision  placement  algorithm. 
Gaussian  noise  was  injected  into  measured  position,  altitude,  and  inertial  sensors.  In  addition  to  sensor  errors,  two 
sources  of  wind  variation  were  added  to  the  simulation;  varying  direction  and  varying  ground  wind  magnitude. 
Wind  was  assumed  to  have  a  constant  magnitude  prior  to  terminal  guidance.  Wind  magnitude  varied  linearly  from 
W  at  the  TIP  altitude  to  WF  at  the  ground.  Prevailing  wind  was  assumed  by  the  system  to  travel  down  range  at  a 
heading  of  zero  degrees  while  the  true  wind  varied  in  its  direction.  For  all  simulations,  the  target  was  set  as  the 
origin.  Sensor  noise  and  wind  variation  statistics  are  listed  in  Table  3.  Dispersion  results  are  shown  in  Fig.  14.  The 
resulting  CEP  is  55.1ft  shown  by  the  circle  and  is  defined  by  the  radius  which  includes  50  percent  of  the  impacts. 
The  results  of  real  drops  made  at  Camp  Roberts,  CA  and  Yuma  Proving  Ground,  Yuma,  AZ  in  the  course  of  2008 
are  discussed  in  Ref. 3 8. 


Table  3.  Error  statistics. 


Parameter 

Mean 

Standard  deviation 

Initial  position  x 

-1500  ft 

250  ft 

Initial  position  y 

Oft 

250  ft 

Initial  position  z 

-3000  ft 

250  ft 

GPS  x  bias 

0.0  ft 

2.0  ft 

GPS  y  bias 

0.0  ft 

2.0  ft 

GPS  x  deviation 

1.5  ft 

0.0  ft 

GPS  y  deviation 

1.5  ft 

0.0  ft 

Altitude  bias 

0.0  ft 

7.5  ft 

Altitude  variation 

1.5  ft 

0.0  ft 

Roll,  pitch  and  yaw  bias 

0.0  deg 

2.0  deg 

Roll,  pitch  and  yaw  deviation 

1.0  deg 

0.0  deg 

u,  v  and  w  bias 

0.0  ft/s 

0.25  ft/s 

u,  v  and  w  deviation 

0.5  ft/s 

0.0  ft/s 

p,  q,  and  r  bias 

0.0  deg 

1.0  deg 

p,  q,  and  r  deviation 

1.0  deg 

0.0  deg 

W 

-10.0  ft/s 

3.3  ft/s 

WF 

-10.0  ft/s 

3.3  ft/s 

Wind  heading  error 

0.0  deg 

15.0  deg 
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Fig.14  Monte  Carlo  simulation  dispersion. 


VIII.  Conclusion 

Using  a  simple  kinematics  model,  an  algorithm  was  developed  allowing  an  aerial  delivery  system  to  predict  a 
when  to  exit  its  loitering  pattern.  The  algorithm  allowed  the  system,  in  the  presence  of  winds,  to  approach  the  target 
at  an  altitude,  so  that  a  final  turn  could  be  completed  with  the  target  being  reached  when  impacting  the  ground. 
Uncertainty  in  winds,  sensors,  and  guidance  resulted  in  the  system  arriving  only  approximately  at  the  predicted  final 
approach  location.  In  order  to  improve  accuracy,  a  final  turn  reference  trajectory  based  on  the  inverse  dynamics  in 
the  virtual  domain  was  generated.  The  method  efficiently  computes  a  reference  trajectory  solution  to  the  two-point 
boundary-value  problem  using  the  current  state  of  the  aerial  delivery  system.  Efficient  computation  allowed  the 
trajectory  to  be  recomputed  during  the  final  turn  adding  robustness  to  changing  winds  and  tracking  errors.  The 
complete  precision  placement  algorithm  including;  trajectory  tracking,  loiter  exit  decision,  and  final  turn,  was 
demonstrated  through  Monte  Carlo  simulations.  The  resulting  algorithm  had  a  circular  error  probable  of  55.1ft  for 
the  example  aerial  delivery  system. 


Appendix:  Model  Parameters 


System  mass  m,  slug : 

Steady  state  aerodynamic  velocity  Va  ,ft/s\ 

0.162 

24.0 

Canopy  reference  area  S  ,ft 2: 

11.0 

Canopy  span  b  ,ft: 

Canopy  chord  c  ,ft: 

Incidence  angle  f  ,  deg.: 

Inertia  matrix  elements,  slug-w2: 

4.45 

2.47 

-12 

4c  =0.3 12,  Iyy 

=0.296,  4  =0.039,  4- =0.022 

Elements  of  the  apparent  mass  matrix  \a  m  ,  slug: 

A  =0.008, 

B  =0.0022, 

C  =0.029 

Elements  of  the  apparent  inertia  matrix  lai  ,  slug-m 2: 

P  =0.04, 

0=0.10, 

R  =0.0018 

x-distance  from  mass  center  to  apparent  mass  center  xBM  ,ft\ 

0.15 

z-distance  from  mass  center  to  apparent  mass  center  zBM  ,ft\ 

-3.64 

Maximum  brake  deflection  d  ,  in: 

5 

Aerodynamic  coefficients:  CD 0  =0.25,  CDa2=  0.12, 

CYp  =-0.23, 

Qo  =0.091, 

Qa  =0.90, 

0,0=0-35,  Cma=~  0.72, 

Cmq=~  1.49, 

Cz^=-0.036,  Clp  =-0.84, 

Cz,  =-0.082, 

ClSa  =-0.0035, 

Cnp  =-0.0015,  Cnp=~ 0.082, 

C„r  =-0.27, 

4^=0.0115. 
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